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We present a theoretical study of inhomogeneous dynamic (resistive) states in a single plaquette 
consisting of three Josephson junctions. These breather states are found to exist in a large range 
of control parameters (dc bias 7, anisotropy r\ and self inductance /3l). We perform an analytical 
analysis of these states, and their stability and resonance behaviour. Numerical calculations con- 
firm our results and show that instabilities arise due to resonant interactions between the breather 
and electromagnetic oscillations (EOs) mediated by the breather state. These instabilities mani- 
fest themselves by various nonlinearities, resonant steps and voltage jumps in the current-voltage 
characteristics. 
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I. INTRODUCTION 

In recent years we witnessed a considerable advanc- 
ing in the area, of study of localized excitations in non- 
linear lattices.trEl These excitations are coined discrete 
breathers and are time-periodic and spatially localized 
solutions of the underlying dynamical equations of mo- 
tion. They generically appear in spatially homogeneous 
systems due to the interplay between discreteness and 
nonlinearity. Discrete breathers were detected experi- 
mentally in different systems, e.g. in weakly couplecLpp- 
tical waveguides^ in the lattice dynamics of solidsjj in 
antiferromagnetx and in systems of interacting Joseph- 
son j unctions .Llili 

Especially the latter example served for many years 
as a well controlled laboratory object to study various 
nonlinear phenomena. Josephson junction systems allow 
for an adequate theoretical description of the underlying 
dynamics. The originally proposed and realized ladder 
geometry of small interacting Josephson junctions with 
an external dc bias became a-sjstem of intense theoreti- 
cal and experimental studies. c3~t3 The dynamic localized 
states in Josephson junction ladders persist due to a well 
known property of a single Josephson junction, namely 
the presence of two different states. These are a super- 
conducting state (zero voltage drop across the junction) 
or a resistive state (nonzero voltage drop across the junc- 
tion). If the junction is underdamped, both states may 
coexist in a wide range of an externally applied dc bias. 
A well known consequence of that is a hysteretic current- 
voltage (I — V) characteristics. 

A breather state in such a ladder is characterized by a 
few junctions being in the resistive state, while the rest 
of all junctions reside in the superconducting one. These 
inhomogeneous dynamic states lead to specific total dc 
voltage drops across the ladder which are used to plot 
I — V characteristics. This method has been successfully 
combined with snapshots made using the low tempera- 
ture scanning laser microscopy technique (LTSJVUr-which 
allows for a spatial resolution of dc voltage drops an Note 



that the experimentally accessible information provides 
only with time averaged data, so the internal dynam- 
ics of breather states is so far not directly observable 
in experiments. The full dynamical picture is however 
more subtle. E.g. the superconducting junctions display 
nonzero ac voltage drops. The amplitude of these ac volt- 
ages decays to zero with increasing spatial distance from 
the resistive junctions. Moreover, a breather state may 
be tuned into resonance with other dynamical modes of 
the system which leads to an increasing complexity of 
the breather dynamics. Experimental studies revealed 
a whole "zoo" of different breather states, instabilities, 
resonances and switching scenarios depending on the pa- 
rameters of the ladder.Lrtll Similar features have been 
found in numerical simulations. 

To clarify the role of control parameters and the cause 
of the primary (generic) instabilities, resonances and 
switchings we carry out a consistent study of a single 
plaquette model containing three Josephson junctions. 
Whereas the occurrence of breathers in an extended lad- 
der may be interpreted as a dynamical breaking of dis- 
crete translational invariance, in the single plaquette 
model the analogy of the translational invariance be- 
comes a permutational invariance. It is this permuta- 
tional invariance which will be violated by dynamical 
states coined breathers in analogy to the ladder case. The 
approach of stu dy i ng inhomogeneous dynamic states in 
reduced systemalall3 has been successfully used in the 
past to get deeper understanding for different breather 
properties as e.g. energy radiation and quantizations 

Such a single plaquette model is the simplest sys- 
tem which supports permutationally broken states in the 
presence of a homogeneous dc bias 7 (Fig.l). It consists 
of two vertical Josephson junctions parallel to the bias 
direction and a horizontal one that is transversal to the 
bias. The Josephson junctions are indicated by crosses. 
In a time-averaged picture we expect only four differ- 
ent states to be realized. The first one is characterized 
by all junctions being in the superconducting state (SS), 
with zero current flowing through the horizontal junc- 
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tion. The second one is the homogeneous whirling state 
(HWS) when both vertical junctions are in the resistive 
state. Due to permutational invariance of this state the 
current flowing through the superconducting horizontal 
junction is still zero. Note that both states are invari- 
ant under permutation, i.e. the exchange of two verti- 
cal junctions. The states of interest are the remaining 
two breather states. Each of these states is character- 
ized by only one vertical junction and in addition the 
horizontal junction being resistive. The second vertical 
junction is in a superconducting state, although there is 
a nonzero current flowing through it. Obviously both 
breather states are not invariant under permutation. 

In the following we will systematically study the depen- 
dence of breather state properties on the control param- 
eters of the system (dc bias, anisotropy, self-inductance). 
Our paper is organized in the following way. In section II 
we present the model, introduce the equations of motion 
and provide with an approximative analytical treatment. 
In section III we introduce a number of numerical meth- 
ods for computing breathers, obtaining their stability and 
following their evolution. With their help we find the do- 
mains of existence, stability and resonances of breather 
states in the space of control parameters and relate them 
to computed I — V characteristics. 



the horizontal and vertical junctions. By making use pi 
the Resistive Shunted Junction model for each junctionO 
we obtain the following set of equations 

Af{<j>\)= 1 + I m , 

Af($0= 7 -/m, (1) 
AT{4> h ) = I m /v , 

The nonlinear operator AT is defined as 

N{<t>) = $ + a<j> + sin(<£) . (2) 

Here, the unit of time is the inverse plasma frequency, 
the dc bias 7 is normalized to the critical current of the 
vertical junctions I cv , and the dimensionless parameter 
a determines the damping in the system. In (^) we in- 
troduce the mesh current i™-that in turn, is determined 
by the flux quantization lawES 

- (3 L l m = <t>\ - 4>l + 4> h , (3) 

where (3l = 2ttI cv L/^q is the normalized self- 
inductance L of the plaquette. Note here, that the volt- 
age drop across a junction is given by <fi. 

B. Breather solutions 




FIG. 1. Sketch of the plaquette with three Josephson junc- 
tions (marked by crosses) in the presence of uniform dc bias. 
Arrows indicate the direction of external current flow (dc bias 
7). 



II. EQUATIONS OF MOTION AND PROPERTIES 
OF BREATHER SOLUTIONS 



A. Equations of motion 

The dynamics of a single plaquette of three Josephson 
junctions (see Fig. I) is determined by the time depen- 



dent Josephson phases of vertical junctions 
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and the 



horizontal junction cfr. We will consider an anisotropic 
system and the parameter of anisotropy is r\ — I c h /I cv , 
where I c h and I cv are respectively the critical currents of 



We are interested in the properties of inhomogeneous 
breather states. In the limit rj = we may exactly 
solve the equations of motion and implicitly obtain the 
breather solution in terms of the resistive state of a single 
vertical junction. In this case the mesh current vanishes 
I m = and the two vertical junctions decouple, while the 
horizontal one follows the dynamics of the vertical junc- 
tions (below we will consider a particular breather case as 
the left vertical junction is in the superconducting state): 



b 1 = arcsm 7 

*a =/(«, 7,*) 
kh 



(4) 



(f> = f(a, 7, t) — arcsm 7 



Here, the function f(a,j,t) describes the evolution of a 
single vertical junction in the resistive state with bias 7 
and. damping a. It is well known that this state exists 
asEI 



> 'V_ = 
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where 7 r is the retrapping current of a single Josephson 
junction. On the other side the superconducting state 
of the first vertical junction (see Eq. (|J)) will persist 
only if 7 < 1. Consequently we obtain that the breather 
solution ([!]) exists for r\ = provided j r < 7 < I . 

For nonzero values of 7/ 7^ exact analytical solutions 
are not available. In the following we obtain approximate 
solutions using a dc analysis, i. e. neglecting the Joseph- 
son phase oscillations. We write the Josephson phases in 
the form: 
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where fl is the breather frequency. Substituting these 
expressions into (|l|) and (||) the mesh current reads 
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Correspondingly the dc bias dependent breather fre- 
quency and the Josephson phase shifts are 
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(8) 
(9) 
(10) 



Here, we may choose c h = without loss of generality. As 
a result we find the I — V characteristics of the Josephson 
junction plaquette in the presence of a breather state: 



7 = Va(l + rj) 



(11) 



The 



where V = <(j> v 2 > = lim^ $ foS(t) - $(0)). 
I — V characteristics of a breather (nil) differs from that 
of the HWS which is 7 = Vol. The breather is char- 
acterized by a smaller voltage drop across the resistive 
vertical junction for a given current value 7. 

As the dc bias current decreases the breather state can 
disappear as a solution of the dynamic equations and 
the system undergoes a transition to the superconduct- 
ing state. By making use of the dc analysis and a stan- 
dard theory of ther*etrapping current in a single small 
Josephson junctionO we expect that the breather state 
disappears as the current flowing through the horizon- 
tal (resistive) junction is equal 7 r . This happens for the 
following value of the external dc bias "fb,r- 



76, r = a(l + 77)- 



(12) 



Remarkably this value of the retrapping current does not 
depend on the self inductance parameter [3^. The depen- 
dence of jb,r on 7/ is shown in Fig. 2. 

For large values of the dc bias the breather state 
switches to the HWS. This switching occurs because the 
left vertical junction can not anymore support the super- 
conducting state. From (^) we obtain that this switch 
should take place at 



I + 77 
l + 2r? 



(13) 
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FIG. 2. Analytic boundaries of the breathers existence in 
(7,77) plane for Pl — 0, a = 0.1. The dc analysis was used. 

The crucial assumption of the above carried out anal- 
ysis is the neglecting of the Josephson phase oscillations. 
Indeed, the inhomogeneous dynamic breather state in- 
duces the electromagnetic oscillations (EOs) in the pla- 
quette and in turn, these EOs change the simple breather 
solution (JsJ) - However, if the breather frequency Q is large 
with respect to the characteristic frequencies of EOs, the 
influence of such a nonresonant interaction between the 
breather and EOs is small. Thus, in the nonresonant 
case the amplitude A of the breather induced EOs is 



V2n 



and correspondingly, it leads to a cor- 
rected prediction of the critical dc bias value "fb,HWS 
where the breather switches to the HWS: 
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lb, HWS 
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The result is shown in Fig. 2. Note that similar existence 
regions for tbs breathers on extended ladders have been 
reported in £3 

The main outcome of the presented analytical treat- 
ment is a triangular shaped region in the space of control 
parameters 7 and 77 where we expect to find breather 
states. Notably the prediction is that for small values 
of damping (a -C 1) the breather should exist for rather 
large values of the anisotropy 77, e.g. for a — 0.05 the 
anisotropy can be as large as 77 ~ 5. Such parameter 
regions are the opposite of what has been initially antic- 
ipated, when the breather states were found by making 
use of the continuation arguments.for the solution (Q) at 
r] = to small nonzero 77 values.tj In this case the hori- 
zontal junctions serve as weak links. Instead we predict 
the existence of breather states for cases of large 77 values 
where the vertical junctions serve as weak links, and the 
horizontal one is a strong link! 
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C. EO frequencies and resonances 

Next we turn to the analysis of the resonant interaction 
between the breather state and small amplitude EOs. 
In order to obtain the characteristic frequencies of EOs 
we linearize the dynamic equations around the breather 
state (|): 

5 v 1 +a5 v 1 +cos(c v 1 )6l = -^-(5 v 1 -5% + 6 h ) , (15) 

PL 

S v 2 +a6 v 2 +co S (nt + 4)S v 2 = ^-(5l-S v 2 +5 h ) , (16) 

PL 

5 h + a6 h + cos(nt)5 h = --^(6 v 1 -5'»+5 h ) . (17) 

Here 8\ 2 (i), 8 (t) correspond to small amplitude EOs ex- 
cited in the presence of the breather state. In the weakly 
damped case (a <C 1) we may neglect the damping term 
and skip the time-periodic terms in ( |l6| ) and (0). As a 
result we derive two characteristic frequencies of EOs: 



(m+-)fi 



(22) 



G=cos(cD^ 

VPL 



(18) 
(19) 

(20) 



These frequencies depend on the anisotropy n and the 
bias current through the phase shift c\ (Eq.(g)). For large 
values of @l the lower frequency w_ ~ } and the large 



frequency u>+ * 
Pl we find oj+ 



y / cos(cf) . In the opposite limit of small 



while 



((r^) 2 -7 2 ) 



2\l/4 



JH -l-MV 

Since the breather is a time-periodic solution and 
the EOs describe small amplitude deviations from the 
breather state, we can expect that the following scenario 
holds. In the case of a proper frequency matching the 
breather may drive EOs into resonance. In turn the in- 
crease in the EO amplitudes will proceed on the expense 
of the breather state. In other words the possibility of a 
breather mediated resonance for EOs turns into a reso- 
nant interaction between the breather and the EOs. In 
such a case we expect in the full dynamical system vari- 
ous instabilities of the breather state and corresponding 
switchings (the voltage jumps) or resonant steps in the 
I —V characteristics to occur. 

Let us classify the possible expected resonance cases. 
The first class of primary resonances may appear if 



(21) 



where m is an integer. These primary resonances corre- 
spond to a matching of the breather frequency or its 
multiples (higher harmonics of the breather state) with 
the EOs frequencies. 

A second class consists of parametric resonances which 
are characterized by 



Finally we may also expect the appearance of combi- 
nation resonances occurring as 



(23) 



In such a case, which is at variance with the first two 
types of resonances, the breather mediates a resonant 
interaction between two different EO frequencies. 

All listed scenarios of resonance involve the breather 
frequency and the EO frequencies. While the approxi- 
mate expression for the breather frequency (|^) is inde- 
pendent on 0l, the EO frequencies depend sensitively on 
the self inductance (|l8|). Consequently the expected po- 
sitions of resonances vary with (3l- Since the existence 
window for breathers (cf. Fig. 2) is expected to not de- 
pend on /3l, certain resonances may simply be shifted 
outside of the breather existence region and thus, be- 
come irrelevant. E.g. for small values of /3l the EO 
frequency w+ becomes very large and we expect then 
only primary and parametric resonances involving uj- . If 
however intermediate values for /3l are chosen, additional 
primary and parametric resonances involving uj + appear. 
But most important we expect combination resonances 
to appear, which can be realized only if the two EO fre- 
quencies, their sum or difference are of the same order as 
the breather frequency. 



III. NUMERICAL ANALYSIS OF BREATHER 
DYNAMICS 

A. Description of numerical methods 

Let us describe the numerical methods of monitoring 
and characterizing breather states. In all cases a Runge- 
Kutta integration method of the differential equations (|l|) 
is chosen. The first way is to fix all parameters including 
the bias 7, to choose some proper initial conditions and 
to integrate the equations for sufficiently long time until 
the system relaxes onto a breather attractor. Once the 
attractor is reached (with sufficient accuracy) all char- 
acteristics of this state may be obtained, including the 
average voltage drop V. After that we may perform a 
small step in 7 and repeat the procedure. Then we find 
the I-V characteristics. This method is especially good 
when approaching an instability, since we obtain explic- 
itly the switching to another state. We also generate 
time-resolved images (movies) of the full dynamical be- 
haviour in order to visually check the realized state or 
monitor the process of switching. The drawback of this 
method is that it is not capable of explaining the nature 
of an observed instability. In addition it may become 
time consuming if the relaxation times become large. 

Another way of studying the system states is to use 
a Newton scheme that allows to iteratively find time- 
periodic states, e.g. breathers. This method is very fast 
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and possesses a high precision. Once a solution is found, 
it may be characterized similarly to the way described 
above. The drawback here is that we can not predict the 
outcome of a switching if a breather state disappears or 
turns unstable. The method is even insensitive to the 
stability property of a state under consideration, so we 
may even miss the switching. 

In order to obtain reliable information about the stabil- 
ity properties of breather states we perform in addition 
a Floquet analysis of small amplitude deviations from 
a breather state which we generated with the help of 
a Newton scheme. We consider the time evolution of 
small perturbations {ej(0), ei(0)} (i = 1,2,3) at t = of 



a breather solution 



^1,2) 1 



l ) over one period Tf, of the 



breather state. The outcome is a map of the space of 
small perturbations onto itself, which is characterized by 
a Floquet matrix M. (of rank 6): 



ei(Tfe) I 1 £i(0) 



(24) 



This matrix is generated numerically. The breather will 
be linearly stable if and only if the matrix M. has no 
eigenvalues with modulus larger than one. The quasi- 
symplectic properties of this map impose the following 
conditions on its eigenvalues: if A is an eigenvalue, then 
A*, e~ aTh /\ and e~ aTh /\* are alsq-eigenvalues (here a* 
denotes complex conjugation of a). El Furthermore there 
is always one eigenvalue At — 1 which corresponds to per- 
turbations tangent to the trajectory of the breather so- 
lution. Consequently there always exists a second eigen- 
value At = e~ aTb . These two eigenvalues will not be of 
further interest. The remaining four eigenvalues can be 
grouped into two pairs each consisting of an eigenvalue 
and its complex conjugated partner. These two groups 
correspond to the two EO frequencies ([II]). In the nor- 
mal case of a stable breather two eigenvalues are simply 
characterized by 

A ± = c - aTb /2^ ± n (25) 

Thus the four relevant eigenvalues will be located on an 
inner circle with radius e~ aTb / 2 < 1 in the complex plane. 
The variation of parameters system leads to a change 
(moving) of Floquet eigenvalues. Resonances correspond 
to collisions of an eigenvalue with its complex conjugated 
partner on the positive real axis (primary resonance), on 
the negative real axis (parametric resonance) or to the 
collision of two different eigenvalues away from the real 
axis (combination resonance) . In all cases colliding eigen- 
values will generally leave the inner circle, one moving 
inside and one outside. Upon further variation of some 
parameter the outer one may cross the unit circle, which 
marks the point of instability of the breather state. For 
small values of a the two events - collision and crossing 
the unit circle - take place close to each other, so that 
a clear connection between instabilities and resonances 
may be drawn. The Floquet analysis thus allows in addi- 
tion to characterize the type of observed instability (res- 
onance) . 



A combination of all three methods - direct numeri- 
cal integration, Newton scheme and Floquet analysis - 
turns out to be a powerful tool which provides with good 
insight of the dynamical processes under consideration. 



B. Domains of breather existence and instability 
windows. 

In order to obtain the breather solution for different 
values of anisotropy 77 and correspondingly, the domains 
of breather existence in the parameter space (7,77), we 
used the following procedure. First, we fixed the value of 
the external dc bias 7 and found a breather solution in the 
form (|j). After that we increased the parameter 77 by a 
small step A77 = 10~ 2 , and the previously found solution 
was used as an initial guess for a Newton method. The 
procedure was repeated for other values of dc bias 7. We 
varied the self inductance parameter (3l while in most 
of our simulations the damping parameter a = 0.1 was 
fixed. 



1. Zero self inductance case ( /3l — 0). 

In the particular case of zero self inductance (5l the 
equations (|l|) are simplified as (j> h — <f>\ — cf>2, and the 
system is characterized by just two degrees of freedom. 
We found that the domains of breather existence in the 
parameter space 7-77 have a triangular form, and it is 
shown in Fig. 3. 

One can also see that the existence domain expands 
in the region of large 77 as the damping parameter a de- 
creases (Fig. 3a and Fig. 3b). The maximum value of 
the anisotropy r) max still allowing to obtain the breather 

state is estimated as r] max ~ — (see inset in Fig. 3b). 

4a 

Note here, as the anisotropy is not large (77 < 1.5 for the 
system with damping a = 0.1) the dc analysis developed 
in the Section II (see, Fig. 2 and Eqs. (12), ( |l3|) ) gives a 
good approximation for numerically simulated domain of 
breather existence (Fig. 3). As the anisotropy increases, 
77 > 1, the region of dc bias 7 where the breather state 
was found, shrinks with respect to the one obtained by 
making use of the dc analysis. 
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FIG. 3. Numerically obtained domains of breather exis- 
tence (thick solid line) and instability windows (dashed) lines 
in parameter space (7,7?). The prediction of the breather 
state instability based on the parametric resonance between 
the breather state and EOs (Eq. ( p^ ) for m = 0) is shown 
by the thin solid line. The self inductance Pl = 0, and the 
damping a is 0.1 (Fig. 3a) and 0.05 (Fig. 3b). The inset 
in Fig. 3b shows the dependence of the maximum value of 
anisotropy r} m ax on the damping a. 

In order to study the interaction of the breather state 
with EOs we carried out direct numerical simulations of 
the system of equations (|l|) and (||) for different values of 
anisotropy. Moreover, we monitored the current-voltage 
{I-V) characteristics of the system decreasing (increas- 
ing) the dc bias by a step of A7 = 0.0005.B The initial 
values of the dc bias 7 corresponded approximately to the 
center of the breather existence domain for the particular 



value of the anisotropy r\. 

We start with the case of a small anisotropy r\ = 0.1 
(Fig. 4). The breather was excited at 7 = 0.5 and it 
persists in a wide range of the dc bias. The breather state 
switches to the HWS as 7 increases beyond some critical 
value, and smoothly transits to the superconducting state 
as the dc bias approaches the retrapping current (Eq. 
(|l2|)). Note here, that for such small values of anisotropy 
the dc analysis (Eq. (|ll|) and dashed line in Fig. 4) 
gives a good approximation practically for a whole I-V 
curve and it means that the resonant interaction of the 
breather with EOs is very weak. 
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FIG. 4. I — V characteristics of a breather state for (3l = 0, 
a — 0.1 and r\ — 0.1. The I-V curve obtained by making use 
of the dc analysis (dashed line) is also presented. The point 
indicates the initial value of dc bias and the arrows shows the 
process of dc bias increase (decrease). 

As we increase the value of anisotropy to 7/ = 1, the 
I-V curve displays peculiar features (Fig. 5). Indeed, 
as the dc bias decreases the I-V curve deviates from the 
linear behaviour ((Eq. (11) and dashed line in Fig. 5), 
displays a cusp at 7 = 0.34 and eventually switches 
to the HWS state. It is interesting that this scenario 
is realized if we start from rather large values of the dc 
bias 7 > 0.33. For lower dc bias values we observe a 
second disconnected branch (the initial value of dc bias 
was 0.27) which displays a back bending and switching to 
the HWS in the I-V curve as the dc bias was increased. 
With decreasing current this breather state transits to 
the superconducting state. 
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FIG. 5. I — V characteristics of a breather state for /3l = 0, 
a — 0.1 and r\ = 1. The dotted line shows the dependence 
of the characteristic frequency combination of EOs on the dc 
bias 7 and the dashed line represent the dc analysis. 

All these features occurring in the current region 
0.28 < 7 < 0.33, can be explained as a parametric 
resonant interaction between the breather state and ex- 
cited EOs. First, we found a good agreement between 
the breather frequency Q = a ci+ri) = ^ anc ^ ^ ne 
frequency of EOs 2w_ = 1.52 ( for 7 = 0.32 ). To 
study this resonant interaction more precisely we carried 
out the Floquet analysis of the breather state. Fig. 6 
illustrates the evolution of the modulus and the argu- 
ment of the eigenvalues of the Floquet matrix M with 
7. The arguments of the eigenvalues can be restricted 
as < Arg(X) < tt. As the dc bias decreases we ob- 
served that for the particular dc bias 7 = 0.34 an eigen- 
value of the Floquet matrix crosses the unit circle at — 1 
(Fig. 6). Thus the original breather state becomes un- 
stable (although it continues to exist as a solution of the 
equations of motion) and this instability is driven by the 
parametric resonance (Fig. 6b) with the EO of the fre- 
quency w_. Note, that as the current decreases further 
this eigenvalue returns into the unit circle, the breather 
state is again stable and thus an instability window of the 
breather is found (the region Ai in Fig.3(a,b)). The po- 
sition of this instability window in the parameter space 
7-77 can be predicted using Eq. (^) or in the form (thin 
solid lines in Fig. 3): 



2V2a 2 (l+r]) 2 
(2m + l) 2 



\ 



(2m + l) 4 



4a 4 (l + 77)2(1 + 2rj) 2 



(26) 
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FIG. 6. a) Modulus and b) argument of the Floquet multi- 
pliers as a function of dc bias 7. The parameters are (3l = 0, 
a — 0.1 and 77 = 1. 

The appearance of parametric resonance is often ac- 
companied by a period doubling of the state mediating 
the resonance, and thus we checked the time dependence 
of the Josephson phase 4>i(t)- The results for two dif- 
ferent dc bias values are presented in Fig. 7. One can 
see that for the dc bias 7 = 0.2965 the period of the 
Josephson phase <p\{t) is two times less than the period 
of <j)\ (t) for the slightly larger value of 7 = 0.2970. Thus, 
the "back bending" (and in same manner the" cusp") of 
the I-V curve (Fig. 5) is related to this period doubling 
process. 
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FIG. 7. Time dependence of the phase <f>\ (t) for two differ- 
ent values of dc bias 7 = 0.2965 (a) and 7 = 0.2970 (b). The 
parameters are /3l = 0, a = 0.1 and rj = 1. 

As the anisotropy rj increases and approaches the tip of 
the domain of breather existence (Fig. 3) the I-V charac- 
teristics start to display different features. First, for such 
a large value of anisotropy (77 > 2) we observed a large 
deviation from the dc analysis (Fig. 8-10). It indicates 
a strong increase of the amplitudes of EOs. Presumably 
this is also the reason why the switching to the HWS oc- 
curs for considerably smaller 7 values as compared with 
the dc analysis (see also Ref . [ |23| ) . This circumstance 
leads to a shrinking of the domain of breather existence 
at the tip. 

We also found that the switching outcome in this re- 
gion of anisotropy upon lowering the dc bias 7 crucially 
depends on 77. Indeed, for rj = 2.1 we observed the 
switching to the superconducting state (Fig. 8) but the 
breather state switches to the HWS (Fig. 9) as the pa- 
rameter r) becomes slightly above the value 2.11. More- 
over, for 77 = 2.3 (Fig. 10) we observed that the switch- 
ing to the HWS occurs via a peculiar intermediate state. 
This state is not periodic anymore (either quasiperiodic 
or chaotic in time), yet it shows up with typical breather 
features when averaged over time. It was stable over the 
whole simulation time of 10 6 . 
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FIG. 8. I — V characteristics of a breather state for (5l = 0, 
a — 0.1 and 77 — 2.1. The dashed line represents the dc 
analysis. 
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FIG. 9. I ~V characteristics of a breather state for /3l = 0, 
a = 0.1 and r\ — 2.11. The dashed line represents the dc 
analysis. 
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FIG. 10. I— V characteristics of a breather state for /3l = 0, 
a = 0.1 and r\ = 2.3. The dashed line represents the dc 
analysis. 

We carried out the Floquet analysis of the breather 
state and found the instability region (region in Fig 
3) that is also due to the parametric resonance between 
the breather state and EOs. At variance with the consid- 
ered above instability window A\, this instability region 
adjoins the lower boundary of the breather existence do- 
main. 



2. Finite self-inductance case ((3l =/= 0) 

The breather properties depend crucially on the self- 
inductance parameter (3^. First, we found that the do- 
main of the breather existence shrinks for small and mod- 
erate values of (3l (see Fig. 11). 




11 

FIG. 11. Numerically simulated domains of breather exis- 
tence in parameter space {~f,rf) for: [3l = (solid line), 
Pl = 0.1 (dashed line), /3l = 0.3 (dotted line) and 
(3l = 0.5 (dot-dashed line). The predicted resonances be- 
tween the breather state and EOs £1 = ui+ for j3 = 0.3 and 
/3 = 0.5 (from top to bottom) are shown by the thin solid 
lines. The damping a is 0.1. 

Moreover, the upper boundary of the breather exis- 
tence domain displays a "saw tooth" form for moderate 
values of the anisotropy (t] > 0.5, see Fig. 11). The 
reason for such a peculiar shape of the breather existence 
domain is the resonant interaction between the breather 
state and EOs (thin solid lines in Fig. 11) leading to 
large amplitudes of EOs and corresponding ac induced 
instabilities of the superconducting state of the left ver- 
tical junction. As the parameter (3l increases to large 
values (see Fig. 12) the resonant interaction between 
the breather state and EOs becomes weak again. In this 
case the domain of breather existence practically coin- 
cides with the analytical predictions in Section II (see 
Fig. 2). 
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FIG. 12. Domain of breather existence in parameter space 
(7, j?) for f3 L = 100 and a = 0.1. 

Also, similar to the case of @l = we found various 
instabilities of the breather state driven by the resonant 
interaction with EOs. These instability windows for a 
particular case of Pl — 1 are shown in Fig. 13. 
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FIG. 13. Domain of breather existence and instability win- 
dows (dashed line) in parameter space (7, rj) for fi L = 1 and 
a = 0.1. The predicted resonances between the breather state 
and EOs are shown by the thin solid lines. The boundary of 
the primary resonance is shown by the dotted line. 
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14. a) I-V characteristics of a breather state for 
/3l = 1, a = 0.1 and rj — 0.5. Thick dashed lines correspond 
to lowering of the dc bias, thin solid lines to an increasing 
of the dc bias. The labels A-,A+,B+ and C relate the reso- 
nances to the observed instability windows in Fig. 13. b) A 
zoom of the I-V curve. The dependencies of EOs character- 
istic frequency combinations on the dc bias 7 are shown by 
dotted lines. The arrows show the various switching outcomes 
as the dc bias increases (decreases) . The time dependence of 
the Josephson phase 0" at 7 = 0.44, V = 2.61 and 7 = 0.25, 
V — 1.482 are shown as insets. 



The breather instabilities manifest themselves by var- 
ious features in the I-V characteristics. For example, at 
the value of anisotropy rj = 0.5 we observed a rich set of 
resonant steps and voltage jumps in the I-V curve (Fig. 



For large values of the dc bias (7 ~ 0.6) a parametric 
resonance of the breather state with high frequency EOs 
fl = 2uj + leads to a resonant step in the I-V curve. More- 
over, the breather state biased to this resonance displays 
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a period doubling bifurcation which manifests through a 
sharp corner in the I-V curve at the bottom of the res- 
onant step (A+ region). As the dc bias decreases to the 
value 7 ~ 0.4 we observe the combination resonance 
51 = LU++LU-. This resonance is characterized by a large 
back bending of the I — V curve and the breather state 
dynamics displays a non-periodic behaviour (see inset in 
Fig. 14 and C region). 

As the dc bias decreases further to the value of 
7 ~ 0.37, a large resonant step in the I-V curve appears. 
It corresponds to the primary resonance with 51 = lo + . 
This resonant breather state switches to the supercon- 
ducting state at the value of 7 = 0.27 as the result of 
an instability driven by the parametric resonance with 
low frequency EOs 51 = 2cj_. These features are clearly 
seen in the dependence of Floquet multipliers on the dc 
bias (see Fig. 15). 
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FIG. 15. a) Modulus and b) argument of the Floquet mul- 
tipliers as a function of dc bias 7. The parameters are pz, = 1 
a = 0.1 and 77 = 0.5. 

The type of resonance can be easily detected from the 
behaviour of the arguments of the Floquet multipliers. 
While unstable solutions (with absolute values of mul- 
tipliers being larger than one) which corresponding to 
parametric resonances are characterized by 'frozen' argu- 
ments at value n, combination resonances are character- 
ized by a merging of arguments over a particular interval 
of dc bias. 

For even smaller values of the dc bias, namely 
0.21 < 7 < 0.26) the breather state becomes again 
stable. For 7 ~ 0.26 this breather state shows up 
with a non-periodic behaviour (see inset in Fig. 14 (b)). 
Most important we observe three new resonances: a com- 
bination resonance 51 = uj + — lu^, a parametric reso- 
nance 51 = |o;_|_ and yet another combination resonance 
251 = lu + + u)_. The breather finally switches to the 
superconducting state at 7 = 0.21. 



IV. CONCLUSION 



We have presented an analytical and numerical study 
of inhomogencous dynamic states (breathers) in a single 
plaquette consisting of three small underdamped Joseph- 
son junctions. Wc found that as the damping a is small 
the breather state as a solution of dynamic equations 
exists in a large region of parameters: dc bias 7 and 
anisotropy r\ (Fig. 2, 3, 12 and 13). For small and mod- 
erate values of anisotropy 77 (77 < 1) the dc analysis 
that completely neglects the breather induced EOs (Eqs. 
(p|)-(|l3|)), results in a good estimation of the boundary 
of breather existence. However, we found that the do- 
main of the breather existence shrinks as the anisotropy 
77 increases, and moreover, for moderate values of the self 
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inductance /3l the upper part of the boundary displays 
a "saw tooth" feature. These effects are due to the exci- 
tation of EOs with large amplitudes. 

The second major result of the presence of EOs is that 
the domain of the breather existence contains a large 
number of instability regions. The physical origin of 
these instabilities is the resonant interaction between the 
breather state and the EOs. It is important that due to 
the presence of two different characteristic frequencies of 
EOs (for the dependence of u>± on the self inductance 
and anisotropy, see Eqs. (|l8|)), a rich set of resonances 
(primary, parametric or combination ones) occurs. Such 
resonant interactions can completely destroy the breather 
state and correspondingly lead to significant shifts of the 
switchings of the breather state to the HWS or the su- 
perconducting state. These switchings appear through 
voltage jumps in the I-V curves. Note here, that in most 
cases we observed the switching to the superconducting 
state as a result of breather instability (parametric or 
combination) but not due to the standard retrapping 
mechanism (Eq. (|l2|)). Similar findings have-been also 
reported for large Josephson junction ladders.E3 

The most interesting observation is that the breather 
state can survive the resonant instabilities by increasing 
its dynamic complexity. This dynamic complexity man- 
ifests itself by resonant steps of various types in the I-V 
curves. Indeed, we observed a resonant type breather as 
a result of a primary resonance with EOs (B + region in 
Fig. 14a), a period doubling process resulting from the 
parametric resonant interaction (Fig. 5, regions A + , A_ 
in Fig. 14) and even a non-periodic (quasi periodic or 
chaotic) behaviour in parameter regions where combina- 
tion resonances appear (C region and the region of small 
current in Fig. 14). 

Finally, all observed effects (switching outcomes, insta- 
bilities) are very sensitive to the damping parameter a, 
and different combinations of the resonances can appear 
as the damping parameter a decreases. 

What can we learn about breather states in extended 
ladder geometries? Our analysis shows that stable 
breathers can be induced for large self inductance. If on 
the other hand one needs switchings between different 
breather states and resonances, intermediate /3l values 
should be chosen. In addition variations of the damping 
(e.g. via temperature control) may increase the complex- 
ity of expected switching scenarios. At the same time it 
should be clear that if already such a small system (with 
three degrees of freedom) shows up with a high complex- 
ity of the phase space structure, the situation for large 
numbers of degrees of freedom as for Josephson ladders 
will be even more complex. 
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